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Abstract: We show how by studying the Earth matter effect on oscillations of solar and 
supernova neutrinos inside the Earth one can in principle reconstruct the electron number 
density profile of the Earth. A direct inversion of the oscillation problem is possible due to 
the existence of a very simple analytic formula for the Earth matter effect on oscillations 
of solar and supernova neutrinos. From the point of view of the Earth tomography, these 
oscillations have a number of advantages over the oscillations of the accelerator or atmo- 
spheric neutrinos, which stem from the fact that solar and supernova neutrinos are coming 
to the Earth as mass eigenstates rather than flavour eigenstates. In particular, this allows 
reconstruction of density profiles even over relatively short neutrino path lengths in the 
Earth, and also of asymmetric profiles. We study the requirements that future experiments 
must meet to achieve a given accuracy of the tomography of the Earth. 
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1. Introduction 

The idea to use neutrinos in order to probe the interior of the Earth was originally put 
forward more than 30 years ago [1] , and since then has undergone a number of improvements 
and modifications [2-15]. In the present paper we offer a new insight into this problem, 
based on some recent theoretical and experimental developments in neutrino physics. 

There are in general two possible approaches to the problem of probing the Earth's 
structure with neutrinos: neutrino absorption tomography [1-6] and neutrino oscillation 
tomography [7-14]. The absorption tomography is based on the attenuation of the flux 
of neutrinos due to their scattering and absorption inside the Earth and so is only sen- 
sitive to the cumulative effect of neutrino absorption and deflection along its trajectory. 
Therefore the absorption tomography cannot give information about the matter density 
distribution along a given neutrino trajectory inside the Earth and requires many baselines 
to reconstruct the Earth's density profile. 

The oscillation tomography of the Earth is possible due to the fact that neutrino os- 
cillations in matter are different from those in vacuum [16, 17]. By studying matter effect 
on neutrino oscillations one can therefore probe the matter density distribution along the 
neutrino path. Being based on an interference phenomenon, the neutrino oscillation tomog- 
raphy has a much richer potential for studying the structure of the Earth. In particular, 
it is in principle possible to use just one baseline and probe the Earth's density at various 
points along the neutrino trajectory. 

In most studies on neutrino oscillation tomography, accelerators were considered as 
the neutrino source. However, solar and supernova neutrinos have a number of advantages 
over the accelerator neutrinos in this respect as the probe of the Earth's interior [13, 14]. 
First and foremost, the Sun and a supernova provide us with a free source of neutrinos. In 
addition, it is very important that, due to the loss of coherence on their way to the Earth, 
solar and supernova neutrinos arrive at the Earth as mass eigenstates rather than flavour 
eigenstates. For oscillations of mass eigenstate neutrinos in a medium, matter effects fully 
develop at much shorter distances than they do for flavour eigenstates [18], therefore by 
using solar or supernova neutrinos one can probe the Earth's density distribution even 
over relatively short neutrino path lengths inside the Earth. Another virtue of studying 
the oscillations of mass eigenstate neutrinos in matter is that in that case asymmetric 
density profiles can be reconstructed [13, 19, 20]. 

Although to first approximation the Earth's density profile can be considered as spher- 
ically symmetric, some deviations from perfect symmetry are possible, especially over rel- 
atively short scales. Exploring such short-scale inhomogeneities would be of particular 
interest from the point of view of the possibility of oil or gas prospecting [12, 13]. As was 
shown in [21], studying asymmetric density profiles through the usual oscillations of flavour 
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eigenstate neutrinos is practically impossible (see also the discussion in section ||). 

Direct inversion of the neutrino oscillations problem in order to reconstruct the matter 
density profile is in general a difficult and subtle mathematical problem - it reduces to 
the reconstruction of the potential of the Schrodinger equation from its solution. In the 
two- flavour case, this in general requires the knowledge of the energy dependence of the ab- 
solute value of one of the components of the neutrino wave function and also of the relative 
phase between the two components [7, 8], which is not measured in neutrino flavour oscil- 
lation experiments. Therefore, one usually has to resort to indirect methods, for example, 
by generating random density distributions and comparing the corresponding predictions 
for oscillation probabilities with simulated data for the "true" profile [10-12]. This is a 
complicated and time consuming procedure of limited accuracy. Indeed, one normally rep- 
resents the matter density profile along the neutrino path as a relatively small number of 
layers of constant and randomly chosen densities, which gives rather poor resolution due 
to the obvious limitation on the number of layers. 

In the present paper we develop a novel direct approach to the neutrino oscillation 
tomography of the Earth. It is based on a recently found simple expression for the Earth 
matter effect on the oscillations of solar and supernova neutrinos in the Earth. The simi- 
larity of this expression to the Fourier transform of the Earth's density profile allows one 
to employ a modified inverse Fourier transformation and reconstruct the density profile in 
a simple and straightforward way. 

The paper is organized as follows. In section [|| we discuss some general features of the 
Earth matter effect on oscillations of solar and supernova neutrinos inside the Earth and 
present the main formulas which will be used for the inversion of the oscillation problem. 
In section |3| we consider the advantages of oscillations of solar and supernova neutrinos 
over the other neutrino oscillations in reconstructing asymmetric matter density profiles. 
In section || we consider the inversion procedure based on the simplest formula for the 
Earth regeneration factor, valid for neutrino path length inside the Earth L <C 1700 km 
(linear regime). We also develop simple iteration procedures which allow one to overcome 
the difficulty related to the lack of knowledge of the regeneration factor in the domain of 
high neutrino energies. In section || we briefly discuss the Earth density reconstruction 
based on the more accurate expression for the Earth regeneration factor, which allows one 
to study longer neutrino path lengths in the Earth (non- linear regime). The requirements 
to the experimental setups which have to be met in order to achieve a given accuracy of the 
reconstructed Earth density profile are considered in section ||. In particular, we discuss 
the effects of finite energy resolution of the detector on the accuracy of the oscillation 
tomography of the Earth. Our main results are summarized and discussed in section [7|. 
Some technical details of our calculations are given in the Appendix. 
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2. Generalities 



In the 3-flavour framework, neutrino oscillations are in general described by two mass 
squared differences, Am|i and Am| l5 three mixing angles, #12, #13 and 623, and the Dirac- 
type CP-violating phase <5cp- For oscillations of solar or supernova neutrinos inside the 
Earth, the third mass eigenstate essentially decouples, and the relevant parameters are 
Am^j, #12 and #13 [22, 23]. The first two of these are determined from the solar neutrino 
data and long-baseline reactor experiment KamLAND [24, 25], while for the mixing angle 
#13 only an upper bound exists. For example, the recent global fit of neutrino data of 
ref. [26] gives for the solar neutrino oscillation parameters the 3a allowed ranges #12 = 
(28.7 -=- 38.1)° and Am^ = (7.1 -=- 8.9) x 1(T 5 eV 2 , with the best-fit values 6 12 = 33.2° 
and Am^ = 7.9 x 10~ 5 eV 2 , while for the mixing angle #13 one finds #13 < 9.1° (13.1)°, or 
sin 6»i3 < 0.16 (0.23) at 90% C.L. (3a). 

Consider a flux of solar or supernova neutrinos arriving at the Earth and traveling 
a distance L inside the Earth before reaching the detector. Due to the loss of coherence 
on their way to the Earth, the incoming neutrinos represent an incoherent sum of fluxes 
of mass-eigenstate neutrinos (see, e.g., ref. [27]). The Earth matter effect on oscillations 
of such neutrinos inside the Earth is fully described by the so-called regeneration factor 
P® e — Pjjg • Here P® e is the probability that a neutrino arriving at the Earth as a mass 
eigenstate v 2 is found at the detector in the v e state after having traveled a distance L inside 
the Earth, and P^J is the projection of the second mass eigenstate onto v e : P^J = \U e2 \ 2 , 
U being the leptonic mixing matrix in vacuum. Note that P^} is in fact the value of P^, 
in the limit of vanishing matter density or zero distance traveled inside the Earth. 

As has been shown in [23], in the 3-flavour framework the Earth regeneration factor 
for solar and supernova neutrinos can be written as 





- cos 2 9 13 sin 2 29 l2 f(5) , 



(2.1) 



where 




(2.2) 



with 



(a;) = \J [cos 2#i2 5 - V(x) /2] 2 + 5 2 sin 2 26 12 
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The effective matter-induced potential of neutrinos V(x) in eqs. ([^) and (|2.3|) is related 
to the charged-current potential Vcc(x) through 



V(x) = cos 2 #13 VccW = c °s 2 #13 

V2G F N e ( 



(2.4) 
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where Gf is the Fermi constant and N e (x) is the electron number density in matter, x 
being the coordinate along the neutrino path in the Earth. The 2-flavour (#13 = 0) version 
of eqs. ([2.1|)-( |2~3| ) was derived in [19], and similar formulas were also found in ref. [20]. 

Equations ( |2.1[ ) and ( ^ ) were obtained under the assumption V(x) <C 26, which is 
very well satisfied for oscillations of solar neutrinos in the Earth. It is also satisfied with 
a good accuracy for supernova neutrinos (except for very high energy ones, which are on 
the tail of the supernova neutrino spectrum). If, in addition, one also requires VL <C 1, 
eq. ( pT2| ) simplifies to [23] 

f(5) = [ L V(y) S m28(L-y)dy. (2.5) 
J 

Equation (^1^) is very suggestive: it has a Fourier integral form and actually means 
that in the small V limit the function f(5) is just the Fourier transform of matter-induced 
neutrino potential V(x) 1 . Therefore, if /(<5) is determined experimentally through eq. (f^H), 
one can employ the inverse Fourier transformation to reconstruct the effective matter- 
induced potential V(x): 

4 f°° 

V( x ) = - I f(5)sm25(L-x)d5. (2.6) 

7T Jo 



The electron number density profile of the Earth N e (x) can then be found from eq. (2.4). 

The Earth density profile could in principle be exactly reconstructed from the solar 
and supernova neutrino data through eq. (|2.6|) under the following conditions: 



1. Eqs. (|2.1|) and (|2.5| ) are exact; 



2. The function f(5) is precisely measured in the whole interval < 5 < 00 (i.e. in the 
infinite interval of neutrino energies < E < 00); 

3. The (5-dependence of the function f(5) is known precisely, i.e the detectors have 
perfect energy resolution, and can determine the energy of incoming neutrinos from 
those of the secondary particles. In addition, the neutrino parameters Am^, #12 and 
#13 are precisely known. 

In reality, none of these conditions is satisfied: eqs. ( |2.1| ) and (^1|) are only valid in 
the limit V(x) <C 25, VL <C 1, the regeneration factor — P^} (and so the function 
f{8)) can only be measured in a finite interval of energies and with some experimental 
errors, the detectors have finite energy resolution and can only give limited information 
on the energy of incoming neutrinos, and the neutrino parameters are only known with 



1 The finite range of integration in (2.5) is related to the fact that the function V(x) is defined on the 
finite interval < x < L. 
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certain experimental uncertainties. In what follows we will study the constraints that 
these limitations put on the accuracy of the reconstructed potential V(x), by relaxing 
conditions (1) - (3) one by one. Conversely, we shall discuss the requirements that are 
put on the experimental installations by the condition of reconstructing the potential V{x) 
with a given accuracy. We will also discuss the ways in which some of the above-mentioned 
limitations can be overcome. 



3. Symmetric versus asymmetric density profiles 

By symmetric density profiles we mean the profiles that are symmetric with respect to the 
midpoint of the neutrino trajectory inside the Earth. They give rise to the potentials that 
have the same property, i.e. 

V(L-x) = V(x). (3.1) 
If the electron number density of the Earth N e was exactly spherically symmetric, the cor- 



responding neutrino potential V(x) would have satisfied eq. (3.1). However, this symmetry 
is only approximate; in particular, it is violated by inhomogeneities of the Earth's den- 
sity distribution on short length scales. Studying these inhomogeneities may be especially 
interesting, e.g., from the point of view of possible oil or gas prospecting. 

Effects of asymmetric density profiles on oscillations of solar neutrinos in the Earth 
have been previously discussed in [13, 19, 20]. Here we give a more detailed discussion of 
neutrino oscillations in asymmetric matter and also compare in this context mass-to-flavour 
and pure flavour neutrino oscillations. 

It is easy to show that in the two-flavour (2f) framework the probabilities of oscil- 
lations between neutrinos of different flavour are the same for the potentials V(x) and 
V(L — x). Indeed, the 2f neutrino oscillation probabilities are invariant under the time 
reversal transformation P a b — > P& a . This follows from the unitarity relations 

Paa Pab — 1 i 

Paa + Pba = 1 , (3-2) 



which enforce P a b = Pb a [21, 28]. On the other hand, for an arbitrary number of flavours, 
time reversal transformation of the probabilities of neutrino oscillations in matter is equiv- 
alent to flipping the sign of the Dirac-type CP-violating phases {<5cp} and replacing the 
potential V{x) with the reverse potential V(L — x) [21]. Since Dirac-type CP-violation is 
absent in the 2f case, from P a b = P^a one immediately finds that 2f oscillation probabilities 
are invariant under the transformation V{x) — > V(L — x). 
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This can also be expressed in the following way. The evolution matrix S of a 2f neutrino 
system is a 2 x 2 unitary matrix which can be written in the flavour eigenstate basis as 

S=( a M (3.3) 
\-/3* a* J 

with \a\ 2 + \P\ 2 = 1. In terms of the elements of S, the oscillation probabilities are given as 
Pab = |5&a| 2 - Hence, time reversal of the evolution matrix (which in the 2f case is equivalent 
to the transformation V(x) — > V[L — x)) reduces to the transposition Sab — ► Sb a , i-e. 

a^a, P^-P*. (3.4) 

Since the transition probabilities P a b = |/3| 2 , as well as the survival probabilities P aa = 
|a| 2 , are invariant under the transformation (3.4), they cannot discriminate between the 
potentials V(x) and V(L — x). This means that flavour oscillations cannot be used for 
a unique reconstruction of asymmetric density profiles: there will always be a two-fold 
ambiguity. 

The situation is drastically different when one considers the transitions between the 
mass and flavour eigenstates, as is the case for oscillations of solar and supernova neutrinos 
inside the Earth. In that case the unitarity conditions read 



Pie + P2e — 1 

p(0) , p(0) 
r le r 2e 



1 , (3.5) 



from which one only finds P2 e — P^e = ~(^~ > ie ~ -^le )> ana - no restrictions on the behaviour 
of the probabilities under the replacement V{x) — * V(L — x) follow. In fact, it was shown 
in [23] that in the 2f case the Earth regeneration factor i-2e — P^^ can be expressed through 
the elements of the matrix S as 

p 2e - pW = cos2#i 2 |/3| 2 + 2sin2^i 2 Re(a*/3) . (3.6) 



While the first term on the right-hand side of eq. ( |3.6j ) is invariant with respect to the 
transformation ( |3.4j ), the second is in general not 2 . Thus, unlike the flavour oscillations, 
the oscillations of mass eigenstate neutrinos into flavour eigenstate ones can be used to 
uniquely reconstruct asymmetric density profiles even in the 2f framework. As was shown 
in refs. [22, 23], in the case of oscillations of solar or supernova neutrinos inside the Earth 
the third neutrino flavour essentially decouples, and to a very high accuracy the problem is 
reduced to an effective 2-flavour one. Therefore, our conclusion about the impossibility of 
using the neutrino flavour oscillations inside the Earth for an unambiguous reconstruction 
of asymmetric density profiles holds also in the 3-flavour case, hence the superiority of mass- 
to-flavour oscillations. This is in accord with the previous findings that matter-induced T 

violation in neutrino flavour oscillations inside the Earth is too small to be measured [21]. 
2 It is only invariant when /3 is pure imaginary, which is the case for symmetric density profiles [21]. 
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4. Linear regime 



The linear regime of the inverse problem of neutrino oscillations in the Earth is based on 
the simple formula ( |2.5| ) for the function /(£), which was derived under the assumptions 

V/2<y<l, VL4Z1. (4.1) 

For the matter density in the upper mantle of the Earth, p ~ 3 g/cm 3 , the second of these 
conditions leads to the upper limit on the allowed neutrino path lengths in the Earth 

L < 1700 km , (4.2) 

which we will now assume to be satisfied. This condition will be relaxed in the discussion 
of the non-linear regime in section |5[ 

4.1 Integration over finite energy intervals 

We shall first assume that eqs. ( |2.1| ) and ( |2.5| ) are exact, but the function f(5) is only 
known in a finite interval [5 m i n ,5 max ] (i-e. in a finite interval of neutrino energies E m i n < 
E < E max ). To study the effects of finite 5 m i n and 5 max , consider the integral of the type 
(|2.6| ) in the finite limits: 

&max 



7T 



/ / (5) sin 2<5(L - x)dS = - / dyV{y) { ^ v -l ^ \ 



(4.3) 

Here we have used eq. ([2.5|), changed the order of integrations and performed the integral 
over 6. 

Ideally, one would like to have S m i n L <C 1 and S max L 3> 1 in order that the integral in 
eq. ( |4.3| ) approach the integral over the infinite interval < 5 < oo in eq. ( |2.6| ) as closely as 
possible. As we shall see, having large enough 5 max in principle does not pose a problem. 
In contrast to this, in most situations of practical interest 5 m i n > L^ 1 , i.e. the condition 
5 m i n L <C 1 is not satisfied, which could be a serious problem. We shall show, however, that 
this difficulty can be readily overcome. 

Let us first study the effect of finite 5 max in eq. (|4.3j). In the limit k — > oo, the function 

sinkx/x goes to 7r5(x) 3 , therefore for 5 max — > oo the upper limit of the integral in ( |4.3| ) 

would yield V(x) — V(2L — x), i.e. essentially the potential V(x) (note that V(2L — x) = 

for all x^ Li). The function = smkx/x for finite > is plotted in fig. [I]. The width 

of its central peak is ~ n/k. With increasing k, the peak becomes higher and narrower, 

and the amplitude of the side oscillations quickly decreases. It is therefore clear that finite 

dmax leads to the finite coordinate resolution Ax ~ Tr/5 max of the reconstructed potential 

V{x) as well as to small oscillations of the reconstructed potential around the true one. 

Note that here S(x) denotes Dirac's delta-function, not to be confused with the parameter 8 defined in 
eq. ((O). 
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Figure 1: Function g(x) = sm(kx)/x 



Large values of 5 correspond to small values of neutrino energy E (see eq. Q2.3] )), 
therefore, the smaller the neutrino energy, the shorter the coordinate scale on which the 
Earth density profile can be probed. This may look somewhat counter-intuitive; however, 
one should remember that we are probing the Earth density distribution with neutrino 
oscillations rather than with direct neutrino-matter interactions. The smaller the neutrino 
energy, the shorter the oscillation length, and so the finer the density structures that can 
be probed. In future low-energy solar neutrino experiments, neutrinos of energies as small 
as a hundred keV to 1 MeV can probably be detected [29-31]; this would correspond to 
Smax — (2 -=-20) x 1CP 11 eV, or Ax ~ TT/S max ~ (3-^30) km, which is a very good coordinate 
resolution 4 . For the neutrino path lengths L > 100 km one finds 5 max L > (10 -j- 100), i.e. 
the condition S max L 3> 1 can be easily met. 

At the same time, as we have already mentioned, having a sufficiently small 5 m i n may 

be a fundamental problem. There are two reasons for that. First, small 5 m i n implies 

large neutrino energies, and there are upper limits to the available neutrino energies. The 

spectrum of solar neutrinos extends up to 15 MeV (hep neutrinos have slightly higher 

energies, but their flux is very small); the average energy of supernova neutrinos is ~ 20 

MeV [32]. Yet, this problem can to some extent be alleviated, at least in principle. It 

might be possible to measure the Earth matter effect for supernova neutrinos of energies 

up to ~ 100 MeV, which are on the high-energy tail of the spectrum, but still not too far 

from the mean energy. Thus, for a nearby supernova and large enough detectors one can 

4 This holds in the idealized case of perfect energy resolution of neutrino detectors. Finite energy reso- 
lution of the detectors is expected to reduce the coordinate resolution of the reconstruction procedure, see 
section 6.2. 
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probably have sufficient statistics. 

The second obstacle is of more fundamental nature. Our expressions ( |2.2| ) and 
are only valid in the approximation V/25 <C 1, which may break down for too small Srnin- 
This gives a lower limit on the values of 5 m i n one can use, depending on the accuracy of the 
reconstructed potential one wants to achieve. In principle, one could attempt at deriving a 
more general expression for f(5), not relying on the approximation V/25 <C 1; in particular, 
it is fairly easy to study the opposite case V/25 3> 1. However, for V > 25 the expression 
for f(5) does not have a simple dependence on the potential V(x), and solving the inverse 
problem of neutrino oscillations becomes a difficult task. 

For an estimate of the constraint on 5 m i n imposed by the condition of small V/25 : 
let us require V/25 m i n < 1/5. For the matter density in the upper mantle of the Earth, 
p ~ 3 g/cm 3 , this gives 5 m i n > 2.8 x 1CT 13 eV, which for Am^ ~ 7.9 x 1CT 5 eV 2 leads to 
Emax i$ 70 MeV. For neutrinos passing through the core of the Earth, the constraints on 
5 m i n and E max will be a factor of 3 - 4 more stringent. 

Let us now estimate the magnitude of 5 m i n L corresponding to L ~ 300 km, which is a 
representative value of neutrino potholing inside the Earth, satisfying condition ( [4.2[ ). For 
solar neutrinos {E max ~ 15 MeV) we find 5 m i n L ~ 2, while for supernova neutrinos, taking 
Emax — 70 MeV, we obtain 5 m i n L ~ 0.43. Thus, except for very small values of L, we have 
to deal with situations when 5 m i n L > 1. 

How large is the error introduced by non-vanishing 5 m i n in the reconstruction of the 
density profile V(x)l To study that, let us consider the integral of the type (|2.6D with 
the finite lower integration limit 5 m i n = L^ 1 . In fig. ^ the lowest curve gives the result 
of such a calculation for the step-function density profile (shown by the dashed line) and 
i5 mai = 300L -1 . One can see several interesting features of the result. First, the deviation 
from the exact profile is relatively small near the detector (x ~ L), but reaches about a 
factor of three far from it (x ~ 0). Second, despite a significant deviation from the exact 
profile V(x), the positions and the magnitudes of the jumps in V(x) are reproduced very 
accurately. Both these features can be easily understood (see sections |4.2| and 4.2.2| below). 



Thus, we have seen that the error in the reconstructed profile due to the lack of the 
knowledge of the Earth matter effect in the domain of low 5 (high energies) can be quite 
substantial. We shall show now how this problem can be cured by invoking simple iteration 
procedures. 

4.2 Iteration procedure I 

Let us study the effect of non-vanishing 5 m i n in more detail. In order to do so, we consider 
the limit 5 max — > oo, which is justified by the preceding discussion. From eq. ([4, 3D one 
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then readily finds 

4 r°° i r L 

V(x) = - f(5)sm28(L-x)d5 + - V(y)F(x, y; 25 min )dy , (4.4) 

7T Jx . VT In 

where the function F(x,y;k) is defined as 

, . smk(x-y) sin A;(2L — x — y) , , . 

F x, y; fc = ^ ^ ^ . 4.5 

x — y 2L — x — y 

It is symmetric with respect to its first two arguments: F(x,y; k) = F(y,x; k). 



Equation (^4) is exact provided that eq. fl2,5|) is exact. By comparing it with eq. (|2.6|), 
we find that the second integral in (^4) can be considered as compensating for an error 
introduced in eq. (|2.6| ) by having a non-zero lower limit in the integral over 5. However, 
this compensating integral cannot be calculated directly because it contains the unknown 
potential V{x). Thus, we have traded one unknown quantity - the function f(5) in the 
domain 5 < 5 m i n - for another. At first sight, this does not do us any good. This is, 
however, incorrect: eq. (|4.4j) allows a simple iterative solution. 



We first note that in the limit <5 min — > the second integral in eq. (4.4) disappears, 
while the first one yields V{x). Therefore, for not too large values of 5 m i n the first term in 
(|4.4j ) is expected to give a reasonable first approximation to V(x). One can then use the 
result in the second integral to obtain the next approximation to V(x), and so on. Thus, 
we define 

4 r°° 

V Q (x) = - f(5)sm25(L-x)d6, (4.6) 



7T 



!U I h 

L 



i r 

I (x) = - V (y)F(x,y;25 rnin )dy , Vi(x) = V (x) + I (x) 

n Jo 



(4.7) 



In-l(x) 



-I V n ^i{y)F{x,y;25 mi/n )dy , V n {x) = V (x) + I n -i{x) . (4.8) 

Jo 

This gives a sequence of potentials Vo(x), V%(x),..., V n (x),... which, for small enough 
bmin, converges to V(x). It should be noted that, while the exact eq. (4^) is independent 
of 

$mim our iteration procedure involves the approximate potentials and so depends on it. 
The smaller the chosen value of 5 m i n , the faster the convergence of V n (x) to V(x); for 8 m i n 
exceeding some critical value (which in general depends on the profile N e (x)) the iteration 
procedure fails. 



5 Its left-hand side is (5 m in-independent, so must be the right-hand side. Actually, by requiring that the 
derivative of the right-hand side of (4.4) with respect to 5 m in vanish, one can recover eq. (2.5) for f(8). 
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This is illustrated in figs. || - 1|. In fig. [| we show the potential V(x) for the step-function 
model of the Earth density profile, along with the zeroth-order reconstructed potential 
Vq(x) and the results of the first, second and fourth iterations V±(x), ^(x) and V±(x) 
iy^{x) is not shown in order to avoid crowding the figure). The calculations were performed 
for 5 m i n = L~ 5 max = 300L . One can see that already the fourth iteration gives an 
excellent agreement with the exact potential. We have checked that for 5 m i n <C L^ 1 already 
the zeroth-approximation potential Vq(x) gives a very good accuracy (see also section 4.2.2| ). 
The wiggliness of the reconstructed potentials in fig. |2|is due to the finiteness of 5 max ; with 
increasing 5 max it decreases. 



•i-J 
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Figure 2: Step-function potential V(x) (dashed line), zeroth order reconstructed potential Vq(x) 
and the results of the first, second and fourth iterations Vi(x), V%{x) and V±(x). The vertical scale 
is that of the corresponding matter density in g/cm 3 , assuming Y e = 0.5 and 9i 3 — 0. The following- 
values of the parameters were chosen: 8 m in = L^ 1 , S ma x = 300L~ . 



If 

fimin exceeds the critical value (which for the chosen profile is approximately equal 
to 2.4L" 1 ), the successive iterations, instead of approaching the true potential, yield the 
potentials which more and more deviate from it. Thus, in this case the iteration procedure 
fails. 

The iteration procedure works very well not only for simple profiles like the step- 
function profile in fig. [2|, but also for more complicated ones. This is demonstrated in 
fig. ||, which is similar to fig. |2|, but was plotted for the PREM-like density profile of the 
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L[km] 

Figure 3: Same as in fig.|], but for PREM-like profile (shown by dashed line), 

Earth [33] . It should be stressed that this profile was used in fig. || for illustrative purposes 
only, since, in the form we used it, it is a realistic Earth density profile for L = 2R§ ~ 12742 
km and not for L = 300 km (for which the Earth density is actually better approximated 
by the step-function profile of fig. |2|). Fig. |] is similar to figs. |2| and || but was produced for 
an asymmetric Earth density profile (shown by the dashed line). It clearly demonstrates 
that asymmetric profiles can also be reconstructed very well, and so the inhomogeneities 
of matter distribution in the Earth can be studied by the method under consideration. 

Finally, we note that if 5 max L 3> 1 but still not very large, one can devise an itera- 
tive procedure correcting for the corresponding error in the reconstructed potential (see 
Appendix |A|). 

4.2.1 Convergence properties of iterations 

We shall now discuss some properties of the iteration procedure ( [4.6[ ) - Q4.8| ). First, we 
note that the function F (x , y; 25 m i n ) is positive definite for all < x,y < L provided that 
5 m in < 2.246L -1 (see Appendix B for the proof). One can then show that for 

5max -> oo and 5 min < min{2.246L~\ 5\} , (4.9) 
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Figure 4: Same as in fig. @, but for asymmetric step- function potential (dashed line). 



where 5\ will be defined shortly, the following string of inequalities is satisfied: 



V (x) < Vx{x) < V 2 {x) < ... < V(x) 



(4.10) 



This can be proven by induction. Consider first eq. ( [4.4] ), From the positivity of V(x) and 
F(x, y; 25 m i n ) it follows that the second integral in ( [4.4] ) (which we denote I(x)) is positive, 
and so the first integral, which is the zeroth-approximation potential Vo(x), satisfies Vq(x) < 
V(x). Generally, the larger 5 m i n , the smaller Vq(x); for large enough 5 m i n the potential 
Vq{x) may even become negative for some values of x in the interval [0, L]. However, if 
S m i n does not exceed certain limiting value 5± (which we assume), the integral Iq(x) defined 



in eq. (4.7) will still be positive. It is this limiting value 5\ that appears in eq. ( [4.9| ). From 
the definition of the integral Io(x) and the obtained condition Vq(x) < V(x) we then find 
< I (x) < I(x). Therefore V\(x) = V (x) + I (x) satisfies V (x) < V\(x) < V(x). Then, 
from the definition of I\{x) we find Iq(x) < I\{x) < I(x), so that V2(x) = Vo(x) + I\(x) 
satisfies V\(x) < V2(x) < V(x). Continuing this procedure, we arrive at eq. ( 4.1 Q|) . 

Thus, under the conditions of eq. (|4.9| ), each iteration produces the potential which 
is larger than that of the previous iteration. The potentials V n (x) approach the exact 
potential V(x) from below and never exceed it. This is well illustrated by figs. [2] - |] (we 
recall that the wiggliness of the curves disappears in the limit 5 max — > oo). Conditions 
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( |4.S| ) are thus sufficient for the convergence of the iteration procedure 6 . 

How can one estimate the critical value of 5 m i n above which the iteration procedure 
would diverge? From the preceding discussion it follows that the convergence of the iter- 
ations to the exact potential relies on the positivity of the "correction terms" I n {x). For 
values 5 m i n > Si the integral Io(x) will become negative for some values of x in the interval 
[0, L], and so for those values of x the potential V\(x) will deviate from V(x) more than 
Vq(x) does. If this propagates into the further iterations, the iteration procedure would 
fail. However, it is possible that even if some Ik(x) are not positive definite, higher-order 
correction integrals I n {x) with n > k are. In this case the iteration procedure would still 
be convergent. 

In practice, it is difficult to determine the critical value of 5 m i n precisely. The closer 
(from below) 5 m i n to its critical value 5 cr a, the larger the number of iterations which 
is necessary to achieve a given accuracy of the reconstructed potential. Therefore, with 
Smin approaching 5 cr u it becomes more and more difficult to check if the procedure would 
converge. In addition, for S m i n close to b cr a the behaviour of the iteration potentials V n (x) 
with increasing n does not depend much on whether 5 ni i n < Scrit or 5 m i n > 5 cr it, so that 
it is difficult to decide if the critical value has already been exceeded. It is therefore 
reasonable to adopt some "practical" definition of 5 cr it, for example, as a value of $min 
for which the number of iterations necessary to reach a 10% accuracy of the reconstructed 
potential reaches 100. For all the density profiles that we studied this value turned out to 
be 5 cr it — 2AL~ l . We elaborate further on the convergence of the iterations in the next 
subsection. 

Let us now return to the question of why the zeroth order potential Vq(x) correctly 
reproduces the positions and magnitudes of the jumps in the exact potential V(x). In 
eq. ( |4.4| ) (which is exact in the linear regime), the second term on the right-hand side is 
the integral over y of the product of the continuous function F(x,y;25 m i n ) and the exact 
potential V(y), which may have discontinuities, but no 5-function type singularities. Hence, 
this integral is a continuous function of x. This immediately means that all possible jumps 
in V{x) are contained in the first integral in eq. (iA), i.e. Vq(x). 



4.2.2 The limit of small 5 m i n L 

It is very instructive to consider the iteration procedure in the limit S m i n L -C 1. Since 
x,y < L, the expression for F(x, y; 25 m i n ) in this case simplifies to 

F(x, y; 26 min ) ~ — <5^ n (L - x){L - y) . (4.11) 



It can be shown that the iteration potentials converge uniformly to V(x) even for larger values of 5, n i n L. 
Indeed, a sufficient condition for the uniform convergence is F 2 (x, y; 25 m i n ) dx dy < ir 2 (see, e.g., [34]), 
which is satisfied for 8 m inL < 2.34. 
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Then from eq. (4JS) we find a very simple result: 



I n (x)~^6 3 



16 

3tt dm '" 1 



V n {y){L-y)dy. 



(4.12) 



It means that the "correction terms" I n (x), which compensate for <5 m j n ^ in the inverse 
Fourier transformation, have a very simple coordinate dependence oc (L — x) for all n and 
scale with 5 m i n L as (5 m i n L) 3 . The same is true for the "exact" correction term I(x) (the 
second integral in eq. (4.4)), which is obtained from eq. (4.12) by replacing in the integrand 
V n (y) by the exact potential V(y). The approximately linear coordinate dependence of the 
deviation of the iteration potentials from the exact one can be seen even for 5 m inL ~ 1 (see 
figs. [|- ||). It explains, in particular, why the deviations of V n (x) from the exact potential 
are small at x ~ L and largest at x = 0. 

Using eqs. ( |4.8| ) and (4.11), it is easy to show by induction that in the limit S m i n L <C 1 

n+r 

/in . . . . \ 

1 



16 



9tt 
3 



u min 2 -' 



/16 , 3 



I 9vr 



L 



L 



V n (x) ~ V(x) - -v 
where the constant vq is defined as 

v = 



16 
9vr 



n+1 



L 



L 



A j\{y){L-y)dy. 



(4.13) 
(4.14) 

(4.15) 



Note that these equations are in accord with eq. fl4.12| ). In the case of matter of constant 
density V(x) = Co = const, one has vo = Co, and eq. ( |4. 14 ) takes an especially simple 
form. 

From eq. (4.14) it follows that for 5 m i n L <C 1 the deviation of the nth-iteration potential 
V n (x) from the exact one scales as [(16/97r) 5^ in L 3 ] n+1 , i.e. the convergence of V n (x) to 
the exact potential is very fast. 

Although eqs. (4.13) and (4.14) were obtained for 5 m i n L <C 1, one can expect that 
they give correct order of magnitude estimates even for 5 m i n L ~ 1. This allows one to 
estimate the number of iterations which is necessary to achieve a given accuracy of the 
reconstructed potential in the whole range 5 m i n < 5 cr u- First, from eq. ( 4.14 ) we find 
that in the case under consideration the critical value of 5 m i n , above which the iteration 
procedure diverges, is 5 cr u = (9-7r/16) 1 / 3 L _1 , so that ( 4. 14| ) can be rewritten as 

V n (x) = V(x) - (3v /2)(5 min /5 CTit ) 3{n+1) (L - x)/L . (4.16) 

Assume now that we want the relative error in the reconstructed potential (which we define 
here as \V(x) — V n (x)\/vo) to be below e. Then from ( 4.16j ) we find that the necessary 
number of iterations is 

ln(2e/3) 



n 



3 ln(S m i n / 5 cr it) 



(4.17) 
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As an example, take e = 0.01. Then for 5 m i n /5 cr it = 0.3 already the zeroth-approximation 



potential Vq(x) has the desired accuracy; for 5 m i n /5 cr it = 0.9 eq. ( 4.17 ) gives n ~ 15, for 
Smin/Scrit = 0.99 it gives n ~ 165, and in the limit 5 m i n — > b CT n one finds n — > oo, as 
expected. It should be remembered, however, that for 5 m i n /Sera ~ 1 eq. ( |4.17j ) gives only 



a rough estimate of the necessary number of iterations because eq. ( [1.16 ) was obtained in 



the limit S m i n L <C 1, which essentially coincides with 5 m i n /5 cr it <C 1. 

We have considered in this subsection the reconstruction of the potential V(x) in the 
limit of small 5 m i n L iteratively only in order to illustrate some general features of the 
iteration procedure. In fact, for 8 m i n L <C 1 it is easy to find the closed-form solution of 
eq. ( |4.4| ) 7 . Indeed, since in this case I (a;) oc (L — x), the potential V(x) can be written as 
V(x) = Vq(x) + C\(L — x) with C± a constant. Substituting this into eq. (|4 . 4|) , one readily 
determines C\, which gives 

V(x) ^ V (x) + C 2 l " 3 " ±j± (4.18) 

with 

° 2 = ihj v ^ L ~y) d y- ( 4 - 19 ) 

Since the function Vo(x) is known, the problem is solved. 
4.3 Iteration procedure II 



In the iteration procedure considered in section |4.2| we assumed that nothing is known 
a priori about the Earth density profile, and the only experimental data available were the 
neutrino data. If some (even very rough) prior knowledge of the matter density distribution 
inside the Earth exists, one can employ a much better and faster iteration procedure to 
reconstruct the Earth density profile. As an example, we consider here the case when 
the average matter density along the neutrino path is known. The corresponding average 
potential 

V = \ f V(x) dx . (4.20) 
L Jo 

is therefore known as well. The new iteration procedure is again described by eqs. (|4.7| ) - 
J4.8D , but the zeroth order approximation potential is now Vq(x) = V = const. In fig. [| we 
plot this potential (horizontal line) and the first iteration potential V\ (x) for the asymmetric 
step- function profile of fig. |] (shown by the dashed line). One can see that, although Vo is 
completely structureless, already the first iteration reproduces the exact profile extremely 
well. This should be compared with the results of iteration scheme I, for which only the 
fourth iteration gives similar accuracy (see fig. |j). 



7 Eq. (i.4) is a linear Fredholm integral equation of the second kind with the symmetric kernel 



F(x,y;2S m i„). In the limit 5 m inL <C 1 the kernel becomes separable (see eq. ( 4.11 )), and the equation 
is trivially solved. 
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Figure 5: Asymmetric step-function profile (dashed line), zeroth approximation potential Vq(x) — 
V and first order iteration potential V\ (x) . Vertical scale and values of parameters are the same as 
in fig. | 



To understand why this iteration scheme works so well, it is instructive to consider it 
in the limit 5 m i n L <C 1. First, we note that for an arbitrary 5 m i n and n > 1 eqs. ( fl.4| ) and 

flU) yield 

V(x) - V n {x) = - [ [V(y) - V n ^(y)]F(x, y; 25 mm )dy . (4.21) 

7T Jo 



In the limit 5 m i n L <C 1, the function F(x, y; 2S m i n ) is given by (4.11); substituting this into 
(4.21), for n = 1 we find 



n.r)-r,(.r)~ (±S 3 min LA^(v -V), 



(4.22) 



where vq was defined in ( [4 . 1 5| ) . For nth iteration (n > 1) we find by induction 

, L — x 



3/16 , , 



L 



-(vo-V). 



(4.23) 



Comparing this with ( |4. 14| ) , we find that in iteration scheme II the difference V(x) — V n (x) 
contains one less power of [(16/97r) £^ iri L 3 ], but is proportional to vo — V rather than 
to vq. We shall show now that for symmetric density profiles (V(x) = V(L — x)) the 
difference vq — V vanishes, and therefore already the first iteration reproduces the exact 
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profile. Indeed, from eqs. ( |4.15| ) and ( [4.20| ) we have 

-L o pL/2 



V 



V = J?\ V(y)[2(L-y)-L)dy = ^ f V(z)zdy, (4.24) 

J J — 



where in the last integral we changed the integration variable to z = L/2—y. For symmetric 
density profiles, V(z) is an even function of z, and the integral in (|4.24| ) vanishes. 

What if the density profile is not symmetric? In that case one can write V(z) = V s (z) + 
V a (z), where V s (z) and V a {z) are the symmetric and antisymmetric parts of V{z). The 
symmetric part V s (z) does not contribute to the integral in (4.24), whereas the contribution 



of V a {z) is small because the inhomogeneities of the Earth density distribution, responsible 
for V a , occur on short coordinate scales 8 . Thus, in scheme II, even for asymmetric density 
profiles the first iteration gives a very good approximation to the exact profile (see fig. |5|) . 

It should be noted that, while the above discussion referred to the case S m i n L <C 1, fig. || 
actually corresponds to 5 m i n L = 1. However, our consideration approximately applies to 
that case as well. For 5 m i n L ~ 1, the integrand of the driving term for the second iteration 
in scheme II will be proportional, instead of V(y)(L/2 — y), to the product of V(y) and the 
function [F(x, y; 25 m i n ) — (1/L) F(x,x';25 m i n )dx']. This function, though not exactly 
antisymmetric with respect to y —* L — y, is almost antisymmetric (for all x £ [0,L]). This 
explains why iteration scheme II works so well even in the case SminL ~ 1. 

5. Non-linear regime 

Our previous discussion of the reconstruction of the Earth density profile was based on the 



simple formula ( |2.5[ ) for the function f(8), which was derived under the assumptions (4.1). 
The second of these conditions led to a rather stringent upper limit on the neutrino path 
lengths in the Earth ( [4.2[ ), or L < 300 km. In order to be able to reconstruct the potential 
V(x) over larger distances, one has to employ an inversion procedure based on the more 
accurate expression ( |2.2| ), which only requires the first condition in ( |4.1| ) for its validity. 
Let us now discuss this improved expression. 

Since the function u)(x) in eq. ( |2.2j ) depends on V(x), we face a non-linear problem now. 
The first condition in ( fd.ip and the non-linearity condition VL > 1 imply 5 m i n L 3> 1. This, 
in turn, means that the matter density profile cannot be found by invoking an iteration 
procedure 9 , and one should resort to different methods of solving eq. (|2.2| ). The simplest 
possibility would be to discretize the problem and reduce the non-linear integral equation 



8 Generally, the contribution of V a to wo — V is of order AVxi/L, where AV is the magnitude of the 
asymmetric inhomogeneity of the potential, and xi is the coordinate scale of the inhomogeneity. While AV 
may be of order of V, Xi/L is always <C 1. 

9 In sec. ^ we showed that in the linear regime the critical values of 8mm, above which the iteration 
approach fails, correspond to 5 m i n L = 0(1). The situation does not improve in the non-linear regime. 
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fl2.2| ) to a set of non-linear algebraic equations, which can then be solved numerically. 
However, eq. (|2.2| ) is a non-linear Fredholm integral equation of the first kind 10 , and 
equations of this type are notoriously difficult to solve. Integral equations of the first 
kind belong to the so-called ill-posed problems: their solutions are very unstable, and to 
arrive at a reliable result one has to invoke special regularization procedures. For linear 



integral equations, such procedures are well developed (see section 6/2 and Appendices C 
and D for more details); however, no universal regularization techniques exist for non-linear 
integral equations of the first kind. The situation is further complicated by the fact that 
for 5 m i n L S> 1 the potential V{x) enters into the integrand of eq. (2.2) being multiplied 
by a fast oscillating function of the coordinate. All this makes the non-linear regime of 
the inverse problem of neutrino oscillations in matter very difficult to explore. This regime 
therefore requires a dedicated study, which goes beyond the scope of the present paper. 



6. Experimental considerations 



Up to now we ignored completely the experimental questions, such as the effects of the 
errors in experimental data and of finite energy resolution of neutrino detectors on the 
accuracy of the reconstructed matter density distributions. We now turn to these issues. 
We shall discuss here only the linear regime of the density profile reconstruction; the 
experimental questions pertaining to the non-linear regime will be considered elsewhere. 

6.1 Effects of experimental errors 

For solar neutrinos, the night-day asymmetry of the signal can be written as [23] 

sin 2 20i2 cos 20i 



A 



ND 



N -D 
'N + D 



'13 



^12 



1 + COS 2012 COS 2012 



(6.1) 



where cos20i2 is cosine of twice the effective 1-2 mixing angle in matter, averaged over the 
neutrino production coordinate inside the Sun [22]. The quantity f(S)/c^ 3 is independent 
of cf 3 in the linear regime (see eqs. ( [2.5| ) and ( |2.4j )). Thus, the relative error of the experi- 
mentally determined value of /(#)/cf 3 is the sum of the relative errors of And, c 2 3 and of 
the 0i2 - dependent expression in the square brackets in eq. (|6T|). The dependence of the 
latter on the error of 0i2 is somewhat involved (mainly, because of the 0i2 dependence of 
the effective mixing angle in matter ^12)- It can be approximated by 



£0i 



1.74 



1.83 [1- (0.4- V av /28f 



(0.4 - V av /28)[1 + 0.4(0.4 - V av /25)\ 



A0 



12 



(6.2) 



We recall that integral equations of the first kind involve an unknown function only under the integration 
sign, whereas in equations of the second kind it is also present outside the integral. 
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where V av is the averaged over the neutrino production region value of the matter-induced 
neutrino potential in the Sun. The values of V av for various components of the solar 
neutrino spectrum can be found in table 1 of ref. [20]. 

The next question is how the error in the experimentally determined quantity /(<5)/cf 3 
affects the accuracy of the reconstructed electron number density profile N e (x). Since, in 
the linear regime, to obtain a solution we invoke an iteration procedure, one might expect 
that the corresponding errors in the reconstructed density profile would accumulate with 
increasing iteration order. This is, however, not the case: at each iteration, the relative 
error of the solution is the same as the relative error of /(5)/cf 3 , and the same applies to 
the exact profile N e (x). This follows from the fact that each iteration profile and the exact 
solution depend linearly on /(<5)/cf 3 , which is a consequence of the linearity of eq. ( f4.4p . 

The main contribution to the error of /(<5)/cf 3 (and thus, of the reconstructed electron 
number density profile) is expected to come from the error in And, which is by far the 
largest one (at the moment, more than 100%). Therefore, an accurate reconstruction of 
the matter density distribution inside the Earth would require very large detectors, capable 
of measuring the solar neutrino day-night effect with an accuracy commensurate with the 
desired accuracy of the matter density reconstruction. 

From eqs. (|2~2| ) and (|2j) it follows that the errors in the parameter Am|i and in the 
energy scale of neutrino detectors go linearly to the shifts in the reconstructed coordinate. 



6.2 Finite energy resolution of detectors 

Let us consider now the effects of finite energy resolution of neutrino detectors. We shall 
be assuming that neutrinos are detected through a charged-current capture reaction, so 
that the energy of the emitted electron in the final state directly gives the energy of the 
incoming neutrino. The electron energy, however, is not exactly measured because of the 
finite energy resolution of the detector, characterized by the resolution function R(T e ,T' e ). 
Here T e and T' e are the observed and true electron kinetic energies, respectively. Because of 
the one-to-one correspondence between the electron and neutrino energies, the resolution 
function can be written directly in terms of the "observed" and true neutrino energies. 
In our discussion it proved more convenient to use 5 = Am|i/4-E instead of the neutrino 
energy E, therefore we will be considering the detector resolution functions expressed in 
terms of 5 and 5' . In many cases the resolution function can be approximated by a Gaussian 

with (5-dependent width, e.g., a(5) = /cq<5 or a (5) = kQyJ8 max 8. 
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When finite energy resolution of detectors is taken into account, the function f(S) in 



the expression for the Earth regeneration factor in eq. (2.1) has to be replaced by 



F{8) = / R(5,5')f(5')d5\ (6.4) 



o 

which describes the smoothing of f(5) with the resolution function. In order to proceed with 
the density profile reconstruction, we have first to extract f(5') from the experimentally 
measured quantity F{8). Once this has been done, one can employ the procedures described 
in sections |2| or ||. Thus, instead of one-step inverse problem, we now have a two-step one. 

6.2.1 Finite energy resolution effects on the Earth regeneration factor 

Let us first discuss the physical effects of the finite energy resolution of detectors on the 
observed Earth matter effect. In ref. [19] it was shown that in the simplified case of the 
box-shaped resolution function with the energy width AE, eq. (|2.2| ) has to be replaced by 



f L ,r, N sinA(5)(L-x) 
Jo A{d)(L-x) 



L 

2 luj{x')dx' 



dx , (6.5) 



where A (5) is the resolution width in terms of 5: A(5) = 5 ■ {AE/E). Eq. ( |6.5| ) differs 
from ( |2.2| ) by the extra factor sin A(5)(L — x)/[A(5)(L — x)] in the integrand. This factor 
equals unity when A(5)(L — x) = (which corresponds to perfect energy resolution or 
x = L), and quickly decreases with increasing A(5)(L - x) (see fig. @). Thus, finite energy 
resolution of detectors leads to an attenuation of the contributions to the Earth matter 
effect coming from the density structures which are far from the detector, the attenuation 
length (distance from the detector) being 

1 E 1 , , E , s 

lau - - sm = -losciE) — . (6.6) 

From eq. Q6.6|) it follows that for the attenuation to be negligible, i.e. for l att to exceed 
considerably L for all energies of interest, one needs 

SmaxL <C . (6.7) 

This can be in conflict with the requirement of good coordinate resolution 5 max L 3> 1, 
unless the energy resolution is extremely good. Thus, finite energy resolution worsens the 
coordinate resolution of the reconstruction procedure. If one completely ignores the fact 
that finite energy resolution of neutrino detectors modifies the observed matter effect and 
simply uses instead of f(S) for the density profile reconstruction, then in the case 

when condition fl6.?1 ) is not satisfied it is essentially the finite energy resolution and not 
the value of 5 max that determines the coordinate resolution of the inversion procedure. 
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Increasing 5 max beyond the limit (6/7) would not lead to any significant improvement of 
the coordinate resolution in that case. 

Eq. (|6.5| ) exhibits a power-law attenuation of the matter density structures far from 
the detector in the case of box-shaped energy resolution function. However, in a more 
realistic case of the Gaussian energy resolution, an even stronger (exponential) attenuation 
results. To show that, we first recall that the function J-(S) is measured in a finite interval 



<5 £ [fiminjfimax]- Although the integral in eq. (6A) extends over the whole interval < 
5' < oo, in fact the contributions to this integral coming from the domains of 5' which are 
far outside the interval [5 m i n , 5 max ] are strongly suppressed because of the presence of the 
resolution function R(5, 5') in the integrand. In particular, if 5 m i n exceeds a few widths of 
the Gaussian resolution function fl6.3| ), one can formally extend the integration in eq. fl6.4| ) 
to negative values of 5' without affecting noticeably the value of the integral. Extending 
the integration to the interval — oo < 5' < oo, from eqs. (|6.4|), ( |6.3| ) and (2.2) one readily 
finds 

L 



F(6) ~ J V(x) e -^ L -*) 2 ° 2 (S) sin 2 J u{x') dx' 



dx . (6.8) 



6.2.2 Undoing the smoothing 

We now turn to the question of how to reconstruct the electron number density profile of 
the Earth from the Earth regeneration factor in the presence of a well-known but finite 
energy resolution. As we already pointed out, if one ignores the fact the observed Earth 
regeneration factor is modified by the finite energy resolution of the detector and naively 
uses J- '(5) instead of f(5) for the density profile reconstruction, the coordinate resolution 
of the reconstruction procedure will in general be very poor. We have checked numerically 
that for 8 max L > 100 and the energy resolution width AE/E = A(5)/5 exceeding 1%, the 
reconstructed profile differs sizeably from the true one. However, if the resolution function 
is well known, one can do a much better job: first, recover the function f(S') from the 
experimentally measured quantity J~(5), and then reconstruct the Earth density profile 
from f(5'). In other words, one could to some extent undo the smoothing of f(5') caused 
by the finite energy resolution of the neutrino detectors and described by eq. ( |6.4| ). This 
is similar to the procedures one has to invoke in various remote sensing problems (e.g., in 
medical tomography). 

In principle, if the resolution function R(5, 5') is exactly known and the function J- {5) 
is precisely determined from the experiment, one could expect that the function f(5') can 



be exactly found from (|6.4D , This is, however, not the case. The point is that eq. (|6.4j) is a 
Fredholm integral equation of the first kind, which belongs to the class of ill-posed problems: 
small variations in J-{S) can lead to very large changes in the reconstructed function f(5'). 
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To illustrate this, let us note that for any integrable R(S, 5') and an arbitrary constant A 



by Riemann-Lebesgue lemma. Therefore, even for very large values of A, changing f(S') — > 
f(S') + As'm(kS') practically does not modify the observable J- {5) if k is large enough. This 
means that small variations in J- {5) correspond to very large high-frequency changes in the 
reconstructed function f(5'). The solution does not depend continuously on the data, i.e. 
is highly unstable. 

The function T (5) is always determined experimentally with some errors; the errors are 
drastically magnified in the process of solving eq. ( |6.4f ) for f(5'). Even if one assumes J~{$) 
to be exactly known, the rounding errors, which are inherent to any numerical calculation, 
would play the same role as the experimental errors and destabilize the solution. Therefore, 
straightforward methods of solving Fredholm integral equations of the first kind (such as 
discretization and reduction to a system of linear algebraic equations) do not work. One 
has to employ some regularization in order to suppress high-frequency noise in the solution. 

There exist many regularization approaches for ill-posed problems and a vast litera- 
ture on the subject. In our calculations we use two popular regularization schemes: the 
truncated singular value decomposition (TSVD) [35] and the Backus-Gilbert (BG) method 
[36-38], which are described in Appendices C and D, respectively. The results of calcu- 
lations for the step-function density profile and the Gaussian detector resolution function 
(|6.3| ) with 10% resolution (ko = 0.1) are presented in figs. |6| and 0. In these figures shown 
are the exact profile (dashed line) and the reconstructed profiles in the case of Gaussian 
resolution: the naive calculation, which used J- (5) instead of f(S'), is shown by the dash- 
dotted curve, while the results of the corresponding regularization approaches are shown 
by the solid curves. In the case of the TSVD regularization, the contributions of the sin- 
gular values Aj < 10~ 10 were truncated (see Appendix C); in the BG approach, 1501-point 
calculation has been used. One can see from the figures that, while the naive calculation 
gives very poor reconstruction of the density profile far from the detector, undoing the 
smoothing caused by the finite energy resolution of the detectors improves the quality of 
the reconstruction drastically. 

Comparing figs. ^ and [7| one can see that the BG method gives in general more accurate 
reconstruction of the profile in the regions that are far from the ends of the neutrino 
trajectory in the Earth, whereas the TSVD method better reproduces the exact profile 
near these endpoints, x ~ and x ~ L, and the jumps of the density (at the jumps, the 
TSVD curve is steeper than that of the BG approach). One can employ both these methods 
of data treating for the same set of experimental data and thus combine the advantages of 
each approach. 




- 24 - 



6 - 



~i 1 1 r 




~i 1 1 r 




100 



200 



300 



L [km] 



Figure 6: Density reconstruction for step-function density profile using TSVD method. The 
dashed line is the exact profile, the dash-dotted curve corresponds to a naive calculation using 3^(6) 
instead of f(5') for 10% Gaussian resolution, the solid curve is the result of TSVD regularization 
for Gaussian resolution, truncation at Xi < 10~ 10 . 



7. Discussion and outlook 



We have developed a novel, direct approach to neutrino tomography of the Earth, based 
on a simple analytic formula for the Earth matter effect on oscillations of solar and su- 
pernova neutrinos inside the Earth. These neutrinos have a number of advantages over 
the accelerator neutrinos from the viewpoint of the tomography of the Earth. Apart from 
coming from a free neutrino source, they are sensitive to the Earth matter effects even if 
the distances they travel in the Earth are as short as ~ 50 - 100 km. In addition, they 
can be used for reconstructing asymmetric density profiles and thus are sensitive to short 
scale inhomogeneities of matter distribution in the Earth, as possibly relevant in prospec- 
tion applications. Neutrino tomography allows the reconstruction of the electron number 
density of the Earth, i.e. is sensitive to its chemical composition. In this respect it can 
nicely complement geoseismological studies, which are based on the measurements of the 
spatial distributions of seismic wave velocities and can only probe the total density of the 
Earth's matter, but not its chemical composition. 

In the linear regime, which is valid for relatively short distances traveled by neutrinos 
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Figure 7: Same as in fig. 0, but for Backus-Gilbert method with N = 1501 points. 



in the Earth (up to several hundred km) , the problem amounts to finding an inverse Fourier 
transform of the Earth regeneration factor. However, to perform this transformation, one 
needs to know the Earth regeneration factor in the whole infinite interval of neutrino 
energies < E < oo, whereas experimentally it can only be measured in a finite range 
Emin < E < E max . Non- vanishing E m i n leads to a finite coordinate resolution Ax ~ 
l sc{Emi n ) of the reconstruction procedure, while finite E max results in a systematic shift 
of the reconstructed profile compared to the true one. 

Future experiments should be able to detect solar neutrinos with energies as small 
as ~ 100 keV, which means that a very good coordinate resolution can in principle be 
achieved; at the same time, finite E max , i.e. the lack of knowledge of the Earth regeneration 
factor in the high energy domain, could potentially be a serious drawback for the neutrino 
tomography of the Earth. We have shown, however, that this ignorance of the Earth matter 
effect at high neutrino energies can be compensated for by making use of simple iterative 
procedures. We have developed two such iteration schemes, one of them not using any prior 
information about the Earth density profile, and the other assuming the average matter 
density along the neutrino path to be known. Both schemes give very good results, the 
second one leading to a faster convergence. 

In order to reconstruct the Earth matter density profiles over the distances of geo- 
physical interest, i.e. exceeding a few hundred km, one needs to solve a complex non-linear 
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problem. It should be noted, however, that even in the linear regime (L < a few hundred 
km) one can obtain very interesting information about the structure of the Earth's crust 
and upper mantle. The maximal depth d of the neutrino trajectory inside the Earth is 



where R® = 6371 km is the average radius of the Earth. For L = 300 km this gives d ~ 1.8 
km, which could be quite sufficient, e.g., for minerals, oil and gas prospecting. Learning 
more about the large-scale structure of the Earth's interior would require calculations in 
the non-linear regime. 

There is a number of very interesting physics and mathematics issues related to the 
problem of neutrino tomography of the Earth, which still are to be explored. We have 
studied in detail the linear regime of the Earth density profile reconstruction; attacking 
the non-linear regime remains a challenge for future investigations. Another issue is recov- 
ering the true Earth regeneration factor, necessary for the inversion procedure, from the 
experimentally measured one, which is smoothed due to the finite energy resolution of the 
detector. Similar procedures are invoked in various remote sensing problems (e.g., in med- 
ical tomography, early vision, inverse problems of geo- and helioseismology, etc.). We have 
used two relatively simple methods of undoing the smoothing effects, the truncated singular 
value decomposition and the Backus-Gilbert method, and only in the linear regime. At the 
same time, there is a variety of other methods developed for solving problems of this kind; 
it would be interesting to study some of them in application to the neutrino tomography 
of the Earth, and also to consider the de-smoothing procedure in the non-linear regime. 

Due to the rotation of the Earth, for a fixed position of the detector solar neutrinos will 
probe the Earth density distribution along a number of chords, ending at the detectors and 
spanning a segment inside the Earth. The size of the segment will depend on the geographic 
location of the detector. Because of the relatively low statistics, only density distributions 
averaged over certain angular bins of incoming neutrino directions will be probed. In order 
to reduce these averaging effects, i.e. to have small angular bins, one would need very large 
detectors and/or large overall detection time. Supernova explosions being one-shot events, 
the Earth tomography with supernova neutrinos will allow the density profile reconstruction 
only along a single neutrino path for each detector and each supernova, depending on the 
detector and supernova positions at the time the supernova neutrinos arrive at the Earth. 

The accuracy of the reconstruction of the Earth matter density distribution will depend 
crucially on the accuracy of the experimental data - of the measured Earth regeneration 
factor and of the neutrino oscillation parameters. It will also depend on the accuracy of 
the reconstruction procedure itself. The largest error is expected to come from the errors 
in the measured Earth matter effect, mainly because this effect is rather small for solar and 
supernova neutrinos. To achieve a given accuracy of the density profile reconstruction, one 




(7.1) 
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should measure the Earth matter effect with similar or better accuracy. The present value 
of the night-day asymmetry of the solar neutrino signal, measured at Super-Kamiokande, 
is 2.1% ± 2.0%(stat) ± 1.3%(syst) [39], i.e. the error constitutes more than 100%. This is 
certainly a very important limiting factor for the reconstruction method described here. 
Future megaton water Cherenkov detectors, such as UNO or Hyper-Kamiokande, are ex- 
pected to reduce this error by a factor of 4 to 10; however, such detectors, unfortunately, 
are not very suitable for neutrino tomography of the Earth. The reason is that they are 
based on the ve scattering process, in which the measured energy of the recoil electron 
gives only loose information on the energy of the incoming neutrino 11 . Detectors based on 
charged-current neutrino capture reactions, such as LENS or MOON [29-31], are probably 
more promising. It is not clear, however, whether sufficiently large (megaton or perhaps 
even tens of megaton scale) detectors of this kind can be constructed, and much work still 
has to be done to assess the feasibility of neutrino tomography of the Earth. In any case, 
studying the Earth's interior with neutrinos is an exciting possibility, which can add yet 
another motivation for very large detectors of solar and supernova neutrinos. 
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A. The case of not very large 8 max L 

If bmaxL 3> 1 but still not very large, one can develop an iterative procedure correcting for 
the corresponding error in the reconstructed potential. 

Let us add to and subtract from eq. ( |4.3| ) the corresponding integral from 5 m i n to <5n, 
where <5o 3> 5 max is a "numerical emulation of infinity" , i.e. a large number for which one 
can neglect the difference between sin(<5ox)/x and ir5(x) (e.g., 5q = 10 3 L _1 ). This gives 

V(x) ~ - / f(6) sin 26 (L - x) d8 + - / V(y){F(x, y; 25 min ) 

7T fx . 7T /n 

+ F{x,y;25 )-F(x,y;25 

max 

)}dy. (A.l) 

If eq. (|2.5| ) is exact, eq. ( |A.l| ) becomes exact in the limit 5q — ► oo. Note that the first term 
on the right hand side of ( |A.l| ) contains integration over 6 from 5 m i n to 6 max , and not from 

11 In ve scattering the neutrino energy could be determined if, in addition to the energy of the recoil 
electron, the direction of its momentum were measured. Precise determinations of these directions are not 
possible with water Cherenkov detectors, but can be realized in ICARUS-type liquid argon detectors [40]. 
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8min to oo. Thus, it properly takes into account that the function f(5) is only measured 
in the energy interval E min < E < E max with E min / 0. 

Eq. ( |A.1| ) can now be solved by iterations, the procedure being quite analogous to the 
one used for solving eq. (4.4) in sections [O or [O. If 5 max L is a very large number, the 



difference between the solutions of eqs. ( [4.4] ) and ( |A.1| ) is numerically insignificant, but it 
can be noticeable when 5 max L is not too large. In our calculation in the present paper we 
always use S max L > 100, so that there is no need to solve the improved equation ( [A. 1| ) . 

B. Properties of F(x,y;28 m i n ) 
Let us prove that 

F(x,y; 25 min ) > for all < x,y < L and 5 min < 2.246ZT 1 . (Bl) 



Indeed, the function F(x,y;k) defined in eq. ( |4.5| ) can be written as F(x,y;k) = g(x' — 
y') — g(x' + y'), where g(z) = sinkz/z, x' = L — x, y' = L — y. For positive k, the 
function g(z) reaches its first minimum at kz ~ 4.493 and is a decreasing function of z 
for < kz < 4.493 (see fig. |]). Since \x' — y'\ < x' + y', F(x,y;k) is positive for all 
< x,y < 4.493//c. Substituting k = 25 m i n and requiring < x,y < L < 4.493/2<5 m j n , one 
arrives at (|B1|). 



C. Truncated singular value decomposition (TSVD) 

Any square-integrable function R(5, 5') can be represented as [35] 

oo 

R(S,S') = X>Ui(<*X*')- (CI) 



i=l 



where Aj are the so-called singular values, which converge to as i — > oo, and Ui(5) and 
Vi{5') are the singular functions. They are orthogonal and can be normalized to satisfy 



(ui-Uj) = I Ui(5)uj(5)d5 = 5ij , (C2) 

and similarly for Vi. There exist standard programs for finding singular values and singular 
functions of square-integrable kernels. 

From eq. d6.4|), which can be written in symbolic form as J- = R- f, and eqs. (|C l|) and 



( |C2|) one finds 

m = ^(ui_Z) Vi{§>) (C3) 



Since Aj — > as i — > oo, the contributions of higher harmonics to f(S') are strongly 
enhanced, which leads to instabilities due to small variations in J- (5). To suppress these 
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instabilities (regularize the solution), one can truncate the series at certain value of i, i.e. 
remove the contributions of very small singular values: 

m „ (C4) 
i=i 

In our calculations we truncate the series when \ become smaller than 1CP 10 . 

A variant of the TSVD method employs a soft truncation, in which in the sum in 
eq. ( |C3|) the following substitution is made: 

1 A,; 



Aj A? + a 2 



(C5) 



where a is the regularization parameter. This procedure is equivalent to using the Tikhonov 
regularization [35]. 

D. Backus-Gilbert method 

This method [36-38] is especially well suited for incorporating the errors of experimental 
data into the inversion procedure. 

We want to solve eq. ( |6.4| ) for the function f(5'). First, we take into account that the 
actual experiments provide us with the binned data, i.e. we will have not a function ^(5), 
but a finite set of discrete values T(5i) = Ti (i = 1, . . . , N). From eq. ( |6.4[ ) we then have 

POO 

Ti = \ R(5 i ,5')f(6')d5' . (Dl) 
J o 

Then, we seek the solution of this equation in the form 

N 

f{5') ~ £>(<T)^, (D2) 

where the coefficients cii(5') are to be determined. This is, actually, the most general form 
for linear inversion. Substituting here T% from eq. ( pi| ), we find 

roa 

f(8')~ 6(6',6")f(6")d6", (D3) 
J( 

where the function 5(5', 5") is given by 

5(5', 5") = ^ ai (5')R(5i,5"). (D4) 



o 



N 



i=l 



From eq. ( p3|) it is seen that in the ideal case the function 5(5' ,5") should coincide with 
Dirac's (^-function. In reality it does not, but we can try to make it as close to the <5-function 
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as possible. The simplest Backus-Gilbert prescription for that is the following. Define the 
spread function r(5') as 



r{5') = / {6' -6"y[5(6',S")YdS" . (D5) 
Jo 

It is non-negative, and would vanish if the function 5(5' , 5") were indeed Dirac's <5-function. 
We now find the coefficients cii(5') by minimizing r(5'), subject to the normalization con- 
straint 

poo 

5(5', 5") d5" = 1. (D6) 



This would lead to 5(5' , 5") being closest to the <5-function for a given choice of the spread 12 . 
The derivation is straightforward, and here we give the results. 
Let us define the matrix S(5') as 



oo 

S ik (5') = I (6' -6") 2 R(5 i ,5")R(6 k ,5")d6" . (D7) 
Jo 

We also define N numbers Ui~. 

/•oo 

R(5 i ,5 1 )d5 1 . (D8) 



If the energy resolution function R(5, 5') is normalized to unit integral, then all Ui in 
eq. ( p§| ) are equal to 1, but in general ui can be different from unity. Minimization of the 
function r(5'), subject to the normalization constraint flD6| ), yields 

N 

"iCO = ir ' (D9) 

J2 uj [5 , - 1 (<5')] jfc 'u fc 

j,k=l 

or, in matrix notation, 

Note that the matrix S is symmetric and positive-definite, so it has an inverse. 



Using eq. (D2), one then finds the approximate solution of eq. (|Dl|): 

N 

f(S') * , (Dll) 

E uj [S , - 1 (<5 / )]jfc«fe 
j,k=i 



or, using the matrix notation, 



2 Other choices of the spread function are also possible. 
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Some comments are in order. To have a good accuracy, one would like the number of 
data bins N to be large, but then S(S') is a very large matrix, and inverting it may be a 
time and memory consuming operation. It is, however, not necessary to invert S(8'). It is 
enough to solve the system of N linear equations 

N 

Ys S *k(S')yk(S') = Ui (D13) 
k=l 



and find an iV-component vector y(S r ), which is a much simpler problem. Then eqs. ( Dll| ) 
and ( pT2|) can be rewritten as 



N 

m - ^ = cow) 

k=l 
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